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ABSTRACT 

Caustics are a generic feature of the nonlinear growth of structure in the dark matter distribu- 
tion. If the dark matter were absolutely cold, its mass density would diverge at caustics, and 
the integrated annihilation probability would also diverge for individual particles participat- 
ing in them. For realistic dark matter candidates, this behaviour is regularised by small but 
non-zero initial thermal velocities. We present a mathematical treatment of evolution from 
Hot, Warm or Cold Dark Matter initial conditions which can be directly implemented in cos- 
mological N-body codes. It allows the identification of caustics and the estimation of their 
annihilation radiation in fully general simulations of structure formation. 
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1 INTRODUCTION 

The idea that the dark matter might consist of a collisionless 
"ga s" of weakly interactin g, neu tral particl es was first p u blishe d 
by ICowsik & McClelland i 19731) and ISzalav & Marxl Jl976t) . 
though earlier discuss i on ca n also be found in the textbook of 
IZeFdovich & Novikivl Jl97lh . These authors proposed neutrinos 
as a promising dark matter candidate, and a claimed measure- 
ment of the electron neutri no mass at very nearly the expected 
value dLvubimov et al.|[l98ll) led to a flurry of interest in neutrino- 
dominated, so-called "Hot Dark Matter" universes. This contin- 
ued until detailed numerical simulations showed the predictions 
for low-redshift s tructure to be quite inconsistent with observation 
dWhite et aljl983h . Attention then shifted towar ds more exotic par- 
ticle candidates for Warm or Cold Dark Matter JPagels & Primackl 
1 1 982l : IPeeblesll 1 982h . 

If a cold collisionless gas evolves from near-uniform initial 
conditions under the influence of gravity, the nonlinear phases 
of growth generically involve caustics analogous to those formed 
when light propagates through a non-uniform medium. This con- 
nection was explored in some depth by R ussian cosmologists 
interested in neutrino-dominated universes jArnoldetal] [l982; 
IZeldovich et~aflll983l) . Caustics are also a very evident feature 
of the similarity solu t ions f or co ld spherical infall p ublished by 
Fillmore & Gold reichl J 1984) an d iBertschingerl (l985). It was an- 
other 15 years, however, before Hoganl teOOlh realised that dark 
matter annihilation could be very substantially enhanced in such 
caustics. He showed that for absolutely cold dark matter the an- 
nihilation probability diverges logarithmically as a particle passes 
through a caustic, and that for realistic dark matter candidates this 
divergence is tamed by the small but finite initial thermal velocities 
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of the particles. He argued that the annihilation radiation from dark 
halos might be dominated by emission from caustics. 

Twenty years earlier lZeldovich & Shandarirj (l 982) had noted 
that thermal velocities limit the densities achievable in dark mat- 
ter caustics, and the first rigor ous calculation of annihi lation rates 
in caustics was carried out by iMohavaee & Shandarirj fcOOrj) for 
Bertschinger's (1985) similarity solution. They found caustics to 
enhance the total annihilation flux substantially in the outer regions 
for plausible values of the initial dark matter velocity dispersion, 
but to be progressively less important at smaller radii. It is unclear 
whether either of these results will apply in general, since the be- 
haviour of the similarity solution is strongly influenced by its spher- 
ical symmetry (which reduces its phase-space dimensionality from 
6 to 2) and by its lack of small-scale structure. 

In the present paper we present a theoretical treatment of the 
growth of structure which shows how the geodesic deviation equa- 
tion can be used to follow local phase-space structure in a La- 
grangian treatment of nonlinear evolution. This formalism is well 
suited for implementation in N-body simulation codes, allowing 
the annihilation signal from caustics to be treated in full general- 
ity provided numerical artifacts from discretisation and integration 
error can be kept under control. We have presented resu l ts from a 
first implementation of this scheme in Vogelsberg er et al.l J2008I) . A 
closely related but some what more complex scheme is described by 
lAlard & Colombil J2005h . Here we complement this work by giv- 
ing a fuller description of the mathematics behind the approach, in 
particular of the regularisation of caustic densities by the finite ve- 
locity dispersion of the dark matter. In the next section we describe 
our idealisation of the initial conditions for structure formation in 
WIMP-dominated cosmologies. Section 3 then presents useful gen- 
eral results for nonlinear evolution from these initial conditions. 
These are used in section 4 to describe the evolution of the local 
phase-space structure, in particular of its projection onto configu- 
ration space, following individual particle trajectories. Caustic pas- 
sages can be identified and the associated annihilation signal can be 
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calculated explicitly. A final section discusses possible future uses 
of this approach. 



2 IDEALISED INITIAL CONDITIONS FOR STRUCTURE 
FORMATION 

In the current standard paradigm for cosmological structure for- 
mation, the dark matter is assumed to be a weakly interacting 
massive particle which decoupled from all other matter and ra- 
diation fields at an early epoch, well before the universe became 
matter-dominated. Since this time, the dark matter has interacted 
with other components only through gravity. In most such mod- 
els there is a period after the transition to matter-domination when 
density fluctuations in all components are linear on all scales, and 
the residual thermal velocities of the dark matter particles are small 
compared to the large-scale velocities induced by density inhomo- 
geneities. In this paper we will take an idealised representation of 
this situation as the initial condition for later evolution of the dark 
matter distribution. 

Thus we assume that at the initial time to the phase-space den- 
sity of dark matter particles can be written as 

f(q,P,t )=p(t )/m P (l+5(g))SS((p-V(g))M, (1) 

where p(t) is the (time-varying) mean mass density of dark matter, 
m p is the dark matter particle mass, q and p are position and ve- 
locity at the initial time and will be used as Lagrangian coordinates 
labelling individual dark matter particles as they follow their trajec- 
tories, 5(g) is the initial linear overdensity at position q, V_(q) is the 
linear peculiar velocity at q and, for growing mode linear perturba- 
tions in the matter-dominated regime, is related to 8(q) through 

V(q) = -Vcj> S(q) cc V 2 0, (2) 

where <f>(q) is the gravitational potential generated by the linear 
fluctuation field, a is the thermal velocity dispersion of the dark 
matter particles at the initial time, and M denotes the standard 
normal distribution in three dimensions. Our assumption that the 
initial density field is linear implies that (S 2 ) <C 1, while our as- 
sumption that the initial dark matter distribution is cold implies that 
<y 2 <§C (|Y] 2 }. Note that the latter condition applies even in the hot 
dark matter cosmology, provided the initial time to is taken suffi- 
ciently late. 

In phase-space the dark matter is thus confined initially very 
close to the three-dimensional "sheet" p — V_(q) and, indeed, its 
distribution becomes exactly three-dimensional in the cold limit 
o~ — > 0. The projection of this sheet onto configuration space (the 
three-density) is very nearly uniform. It proves useful to work in 
initial coordinates where the velocities at each point are taken rela- 
tive to V_(q). We therefore define 

A(q,p)=p-V(q), (3) 

and we approximate the initial phase-space density of the dark mat- 
ter as 

f(q,A,t )=f Af(A/a), (4) 

where we neglect the spatial modulation by the factor (1+8) so that 
/o = p(to)/m p and the phase-space density becomes independent 
of q. Because V_ is the gradient of a scalar field the second-rank 
tensor dA/dq is symmetric and the transformation o f variables 
(q,p) «-» (q,A) is canonical ( IBinnev & Tre maine 2008]). This will 
important below. 



3 EVOLUTION OF THE DARK MATTER 
DISTRIBUTION 

We will assume that dark matter particles interact purely gravita- 
tionally at all times of interest. Each then moves independently of 
the others subject only to the collective gravitational potential. Par- 
ticle accelerations depend on position and time, but not on velocity, 
and trajectories can be derived from th e time-dependent Hamilto- 
nian of the system (e.g. Peebles 1980). The evolution of the dark 
matter distribution is thus a Hamiltonian flow. A thorough develop- 
ment of the properties of su ch flows can be found in Appendix D of 
iBinnev & Tremainel d2008h . We characterise the phase-space posi- 
tion of a particle at time t by its position x and its peculiar velocity 
v. In a Hamiltonian flow, particle trajectories through phase-space 
never intersect, so we can write the phase-space position at time t as 
a unique and invertible function of the initial phase-space position, 
i.e. 

x=x(q,p,t), v = v(q,p,t) 
with the well-defined inverse relation 

q = q{x,v,t), p=p(x,v,t). 

This is a standard Hamiltonian map so the transformation (q,p) <-> 
(x, v) which it defines is canonical, incompressible and symplectic 
(see IBinnev & Tremainel [20081 . for the mathematical significance 
of these terms). As noted above, the transformation of variables 
(q,p) <-> (q,A) is canonical, so we can restate the (canonical) 
relation between initial and final configurations as 

x = x(q,A,t), v = v(q,A,t) 

with 

q = q(x,v_,t), A = A(x,v,t). 

This simplifies the description of our problem considerably. 

The collisionless Boltzmann equation is an immediate conse- 
quence of the incompressibility (in 6-D) of Hamiltonian flows - 
phase-space density is conserved along every trajectory in the flow. 
Eq. (O then provides a complete formal solution for the evolution 
of the phase-space density distribution of the dark matter distribu- 
tion: 

f{x,v,t)= f M{A{x,v,t)/a). (5) 

The maximum phase-space density at time t is fo, and this den- 
sity is achieved everywhere on a 3-dimensional subspace defined 
implicitly by A(x, v_, t) = 0. Phase-space densities are only signif- 
icantly different from zero at points which are sufficiently close to 
this subspace which we refer to below as the "central sheet" of the 
phase-space distribution. The geometry of this sheet is very simple 
at early times: its projection onto 3-space is (approximately) uni- 
form and only one point of the sheet projects onto each x-position. 
Non-linear evolution stretches and folds the sheet, but does not tear 
it. It can then pass through a given ^-position multiple times, pro- 
ducing a series of streams, each with a different velocity v. Caus- 
tics arise on the boundaries between regions with different num- 
bers of streams. Figure 1 illustrates this situation for an analogous 
1 -dimensional system. 

It is important to realise that the solution of Eq. $5$ depends 
only on the assumed initial condition and on the fact that the dark 
matter obeys the collisionless Boltzmann equation. Thus it holds in 
the absence of any symmetry and during strongly non-equilibrium 
phases of evolution. In addition, it does not assume the gravitational 
potential to be generated by the dark matter alone, so it is valid, for 
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Figure 1. Illustration of our calculations for an analogous 1 -dimensional 
system. At the initial time to , the dark matter phase space density is non- 
zero only in regions of phase space close to the central sheet p = V(q), in- 
dicated by the solid curve in the upper left diagram. The shaded regions sur- 
rounding this curve indicate the 1 and 1-a regions of the (Gaussian) phase 
space density distribution. For convenience, we change the initial velocity 
coordinate in phase space to A = p — V(q). As shown in the upper right di- 
agram, the equidensity contours of the phase space density then correspond 
to A = const., and the maximum phase space density occurs on the A = 
axis. Dynamical evolution distorts these initial phase-space distributions ac- 
cording to the Hamiltonian flow (q,p) <-» (x, v), producing a phase space 
distribution at a later time t which is indicated schematically in the lower 
diagram. The collisionless Boltzmann equation guarantees that the phase 
space density is preserved by this map, as indicated by the shaded regions. 
Our plots indicate how two points R and C, corresponding to two different 
dark matter particles, are transformed by these maps. The (1 -dimensional) 
space density in the neighborhood of each is given by projecting the phase 
space density down the velocity axis, and is uniform in the initial space. At 
time t, R is at a regular point of its trajectory (dA/dv ^ 0), whereas C 
is passing a caustic (dA/dv = 0). Clearly the space density at R depends 
on the local slope of its A = const, line, while at C it depends on the 
curvature of this line and on the offset of C from the central sheet of its 
stream. 



example, in the inner regions of galaxies, where the gravitational 
effects of the baryonic components are dominant. These influence 
the trajectories of individual dark matter particles, and so the details 
of A(x, v, t), but they do not affect the Hamiltonian nature of the 
flow or the validity of Eq. {3J. 

The Hamiltonian nature of the flow has a number of con- 
sequences for the map {q,A) <-> (x,y_). If we define 6- vectors 
w = (q, A) and W = (x, v), then the 6-tensors 



= ^dW 

aw 



and W=%!L 
dW 



satisfy the relations 



D D' = D> D = I, det(L>) = det(D') = 1. 



(6) 



The first relation merely states that the backward transformation re- 
verses the forward transformation, so the matrix corresponding to 
the former is the inverse of that corresponding to the latter. The sec- 
ond relation states that both matrices have unit determinant so that 
the transformations are volume- and orientation-preserving. The 



conservation of phase-space density expressed by the collisionless 
Boltzmann equation is a consequence of this second property. Fur- 
ther important properties follow from the fact that these matrices 
are symplectic. In particular, the velocity and space parts of the 
forward and backward transformations are related by 

dg _ dA % _ dg 82 _ dA dv _ dq 
dq dv ' dA dv ' dq dx ' dA dx 



(7) 



where in each equation the partial derivative on the l.h.s. refers to 
the forward map so that the independent variables are q and A, 
while the partial derivative on the r.h.s. refers to the reverse map 
so that the independent variables are x and y_. We will use some of 
these relations later. 



4 VARIATION OF THE 3-DENSITY ALONG PARTICLE 
TRAJECTORIES 

The differential annihilation probability for an individual dark mat- 
ter particle depends on the (velocity-dependent) annihilation cross- 
section cr x (v) and the local phase-space distribution as 



-^jfev,*) = j ' d 3 v' f(x,v_,t)cr x (\v_-v\) \v_ -v\. 



(8) 



In the following we will assume that a x (v) oc 1/v as is the case 
for many (but not all) dark matter candidates. This relation then 
simplifies to 

^(x,v,t) = (a x v) J d\'f(x,vL,t) = p(x,t), (9) 

where p(x, t) is the local 3-space mass density of dark matter and 
{o- x v} is the thermally averaged velocity-weighted annihilation 
cross-section. The total annihilation probability over a finite time 
interval is thus simply obtained by integrating the local dark matter 
3-density along the particle trajectory. This density is made up of 
two distinct components, one due to particles which are part of the 
same stream as the particle in question, and so were its neighbours 
in the initial conditions, and one due to particles in other streams, 
which typically originated in distant parts of phase-space. Caustics 
arise in the first of these two components and so we will concentrate 
on it in the following. 

For our assumed initial condition, each particle can be speci- 
fied by its initial phase-space position (q ,A p ), where \A p \ ~ a 
is very small. The subscript p here identifies that the coordinates 
belong to the specific particle under consideration; it has nothing to 
do with the initial phase-space coordinate p. The particle's later tra- 
jectory is then x p (t ) = x(q ,A p ,t),v p (t) = v(q p ,A p ,t), and the 
3-space stream density at its position can be obtained by integrating 
the phase-space density over all velocities: 



i \A (x„, v,t) I 2 
•- exp ( - '-^t'-' n ] . (10) 



2a 2 



The velocity integral is restricted to velocities such that q (x_ pl y_, t) 
remains in the neighborhood of q . Since a is very small, we can 

simplify by carrying out a Taylor expansion of A (x, v, tj around 



(11) 



where Sv — v — v p and the partial derivatives are evaluated at 

(x p ,v p ). 
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4.1 Densities at regular points 

At almost all points of the particle's trajectory the linear map in 
the second term on the r.h.s. of Eq. Jl lb is non-singular (i.e. the 
determinant of the corresponding matrix is non-zero). These will 
be called "regular points" of the trajectory in the following. At such 
points there exists a small value of 8v, say 5v c , for which the sum of 
the first two terms vanishes; (a; , v + Sv c ) is then the intersection 
at time t of x = x p with the central sheet of the stream to which 
our particle belongs. The integral for p s (x p ) becomes particularly 
simple if we centre the velocity integration on this point. To lowest 
order, we have 



\A\ 



,dA fdA 

~ dv \ dv 



5v' 



(12) 



where Sv' = v — y_ p — 5v c . The tensor product (dA/dv)(dA/dv) T 
is clearly symmetric and must have positive eigen-values. Without 
loss of generality, let us define the principal axes in velocity space 
so that these eigen-values are s 2 ^ s 2 ^ s§ > 0. The condition 
that our point be regular forces the final strict inequality, since an 
obvious consequence of Eq. d 1 2b is 



det -= 



<:)r 



sis 2 s 3 . 



(13) 



In this velocity frame the integral in Eq. d 1 Ob takes a simple form 
which we can easily evaluate, 



Ps(x p (t)) 



fom p 

(2-k 0-2)3/2 



1 



fom p 



S1S2S3 



= fom p 



det -= 



di< 



(14) 



= fomp 



aq 



The last equality follows from one of the relations listed in Eq. {7J 
and demonstrates explicitly that the stream density obtained here is 
identical to that obt ained by forward integr ation of the geodesic de- 
viation equation in lHelmi & White! i 19991) and jyogelsbergeretal] 
(2008). This is a manifestation of the fact that the local velocity dis- 
tribution in a stream is distorted in a way which exactly mirrors that 
of the density field. To lowest order, the map rotates each infinites- 
imal Lagrangian volume and then stretches it by different amounts 
si, S2 and S3 along three orthogonal axes. (Note that the Si can be 
negative.) The (initially isotropic) velocity distribution at the cen- 
tral point is compressed by these same factors along the same set 
of axes. 



4.2 Densities at caustic crossing 

As a particle moves along its trajectory, it will occasionally pass 
through discrete points where the determinant of OA/ dv vanishes 
and the corresponding linear map becomes singular. At such points 
at least one of the stretch factors must vanish, and Eq. dl4b predicts 
an infinite stream density. These points are the caustics of the map. 
In the following we will neglect higher order singularities, where 
two or more stretch factors vanis h, and assume s 2 ^ s| > 0, s| = 
at such caustic crossings (see lTremaind[l999l. for discussion of 
higher order singularities). 

To simplify the integral in Eq. J 10b it is useful to take coordi- 
nates in ?;-space along the principal axes of (dA/dv)(dA/dv) T . 
Unit vectors along the axes corresponding to s 2 and s\ are ro- 
tated into a pair of orthogonal vectors in yl-space by the linear map 



dA/dv, and in addition are stretched by factors si and S2. We use 
these rotated directions to define our 1 and 2 axes in A-space. All v- 
vectors are then projected onto the 1-2 plane in A-space. To lowest 
order, Eq. dl lb becomes 

A(x p , v, t) — si 8v[ e x + S2 5v' 2 e 2 



+ 



A p ,3 + S3 5v 3 + i g V3 2 J g 3) (15) 



where the e.^ are unit vectors along the coordinate directions in A- 
space. We have shifted the origin in ^-space in order to simplify 
the coefficients of e t and e 2 . By allowing the shift to depend on 
Svs, all terms independent of 8v\ and Sv2 can be removed; remain- 
ing second-order terms are then small compared to the linear term 
which is retained. A similar manipulation is not possible for the co- 
efficient of e 3 because its linear term vanishes at caustic crossing 
(i.e. S3 = when t = t c ) and we must retain both constant and 
quadratic terms. Note that the only quadratic term we need to re- 
tain is that involving Sv$ alone, since all others are subdominant to 
the linear terms involving 5v[ and Sv' 2 or are removed by the ori- 
gin shift made for each value of 8V3 . Inserting this expression into 
Eq. J lOt and setting S3 = (i.e. t — t c ), the integrals over vi and 
V2 can be carried out as before and we are left with 

Ps{X p {t c )) = 



\s\si\ (2ttct 2 ) 1 /2 

+ 00 

J dvz exp 

-00 

fomp T' (±\A p , 3 \/ff) 



1 

2^2 



sis 2 



A P ,3 + - -5-5- OV3 



cxp(-^ 3 /20- 2 ) (16) 



icrdMa/c^Ij 1 ^ 

where the dimensionless function V' is of order unity and is defined 
by 



r'PO = 



(2tt) 1 /2 



dy exp (X 2 /2 - (X + y 2 /2) 2 /2) . (17) 



The sign of the argument of T' in Eq. d 1 6b is positive when A Pt s 
and d 2 As/dv3 have the same sign and negative otherwise. Note 
that r' can be expressed in terms of a modified Bessel function, but 
we forgo the details here. 

As before, the mass density at the particle's position is re- 
lated to the initial density fom p by the product of three dilu- 
tion/compression factors. The factors associated with the two axes 
in the plane of the caustic correspond to those we obtained above 
for regular points of the trajectory; caustic formation does not sig- 
nificantly effect the distortion of the dark matter distribution in di- 
rections parallel to the caustic. The compression/dilution in the di- 
rection perpendicular to the caustic is of different form, depending 
explicitly on how cold the initial conditions were (i.e. on the value 
of a) and on how close the particle is to the central sheet of its 
stream (i.e. on the value of A p ^/a), as well as on the overall struc- 
ture of the Hamiltonian flow, as encapsulated by the second deriva- 
tive, d 2 Az/dvi, (x p , v p , t c ). Note that as the initial conditions are 
made colder, the maximum density achieved during caustic passage 
increases as <t -1 ^ 2 . 

For particles which lie away from the central sheet of the 
stream, A p ,s is non-zero. If A Pi 3 and d 2 As/dv 2 have the same 
sign, the maximum of ps(x p (t)) then occurs either slightly be- 
fore or slightly after t c . This is because a small but non-zero value 
of S3 allows the linear term in the coefficient of e 3 in Eq. dl5b 
to be significant. For one particular S3 the resulting quadratic has 
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an extremum at A3 — 0. This value of S3 obtains at the moment 
when our particle is spatially coincident with the caustic of the cen- 
tral sheet of its stream, and so sees the corresponding density (see 
Fig. 1). The exact value of the maximum density is not important 
in practice, since it appears only in the argument of a logarithm 
when we estimate the total annihilation radiation from caustic pas- 
sages. For simplicity, we will set A p .s = when we use Eq. d!6t 
below to estimate the maximum density during a caustic passage. 
In this case, V' can be expressed in terms of the standard complete 
T-function, r'(0) = 2 5/4 r(5/4)/0F. 

4.3 Integrating annihilation rates through caustics 

The formalism developed above is particularly powerful when em- 
bedde d in a high-resolu t ion si mulation of cosmic structure forma- 
tion. IVogelsberger et alJj2008h showed how the geodesic deviation 
equation can be integrated in parallel with the TV-body equations 
of motion to give the full phase-space distortion tensor D along 
the trajectory of each simulation particle. This is then sufficient 
(through Eq. Q) to give all the first order derivatives of the for- 
ward and backward Hamiltonian maps we have been discussing. In 
particular, this allows the local stream density to be calculated at 
all regular points of each trajectory using Eq. d!4t . and this can be 
inserted into Eq. ([9} to obtain the intra-stream annihilation rate at 
such points. 

Caustic crossings can be recognised in an A-body integration 
by the change in sign of the determinant of dx/dq. Let us denote 
by (si, S2, S3) and (s^, s 2 , s 3 ) the stretch factors at the beginning 
and end of a timestep during which a caustic crossing is detected. 
Then si ~ s[ and S2 ~ s' 2 , while S3 and s 3 should be much 
smaller in absolute magnitude and should have opposite sign. A 
good approximation to the evolution of the distortion tensor during 
the timestep is then obtained by assigning mean values to si and s 2 
and assuming S3 to vary linearly between its endpoint values. Ac- 
cording to the arguments of §4.2, the maximum value of p s during 
the timestep is then well approximated as 

2 5 / 4 T(5/4) f m p 



pmax 



(18) 



Away from caustic crossing (but still within the timestep) Eq. d!4t 
can be used to obtain p s which then varies inversely with S3 and 
so with \t — t c \, the time to caustic crossing. This suggests the 
following approximation to the variation of stream density within 
the timestep: 

f m p At 



PsilL v {t)) 



\SlS 2 \ I S3 -S 3 | 



((* 



t c ) 2 +T 2 y 1/2 - 



(19) 



where At is the length of the timestep and T is chosen so that 
Ps(x p (t c )) = /5 max . The shape of this function in the neighborhood 
of its maximum is chosen purely for convenience, but the wings and 
the maximum value itself are correct. If we integrate Eq. ((9} over 
the timestep using this formula we obtain 



AP = 



(crxv)fo At ^ f 2 g / 2 r(5/4) 2 js3S 3 | 



sis 2 s 3 



\d 2 A 3 /dv : 



(20) 



where we have used the fact that T <C At and, for simplicity, we 
have assumed that the caustic occurs well away from either end 
of the timestep. Since a single simulation particle represents many 
dark matter particles, AP must be multiplied by m a i m /m p (where 
m s i m is the simulation particle mass) to obtain the total number of 
annihilation events. 

The properties of the dark matter appear in Eq. d20t through 



the annihilation cross-section, and through the parameter a which 
appear s in the argum ent of the logarithm. Thus we recover the re- 
sult of Hogan (2001) that in the cold limit the annihilation lumi- 
nosity of a caustic is logarithmically divergent. For given parti- 
cle physics parameters, an integration of the A-body and geodesic 
deviation equations provides all the quantities needed to calculate 
AP for each caustic passage, with the important exception of the 
second derivative d 2 A3/ dv 2 . Since this quantity appears only in 
the argument of the logarithm, it is sufficient to estimate it to or- 
der of magnitude, provided the estimator chosen has no strong bias 
when averaged over many caustic passages. 

The condition for caustic passage, det(dA/ dv) = places 
no constraint on the values of the second derivatives, so we can es- 
timate the size of a typical component of d 2 Aj dv 2 from the avail- 
able Galilean-invariant quantities. These are the components of the 
6-dimensional distortion tensor D and the particle acceleration a. 
For example, we can note that a(d 2 A/dv 2 ) (dx/dA) is dimension- 
less and is expected to be of order unity, so that the size of our de- 
sired component of the second derivative can be estimated as the 
inverse of the product of the r.m.s. sizes of the components of a 
and of the components of dx/ OA. This assumption can be checked 
in simple 1 -dimensional cases, and we will give an example based 
on the similarity solution for spherical collisionless collapse in a 
future paper. In this case, at least, the estimate we advocate here 
works remarkably well. 



5 DISCUSSION 

In this paper we have developed the mathematical background to 
enable a relatively precise evaluation of the annihilation radiation 
from dark matter caustics in fully general simulations of the non- 
linear growth of structure. Our scheme allows the annihilation rate 
to be integrated along the trajectory of each simulation particle, in- 
cluding correctly the contributions from all the caustics in which it 
participates. Typically each particle experiences several such caus- 
tic passage s in each orbit around the dark matter halo in which it 
resides (see Vogelsber ger et alj2008h . In order to include correctly 
the annihilation rate between particles which are members of differ- 
ent streams, it is necessary to estimate a local coarse-grained den- 
sity at the position of each particle, and to add in the contribution 
due to streams other than its own. This can be done, for example, 
using the SPH technique, since the smoothing this introduces does 
not bias the luminosity predicted from inter-stream annihilations. 
When a particle passes through a caustic occurring in a stream other 
than its own, the time-integrated annihilation probability is still cor- 
rectly reproduced in the smoothed system. 

By implementing these methods in high-resolution simula- 
tions of galaxy formation it should be possible to achieve a com- 
plete numerical description of the expected annihilation radiation, 
limited only by the ability to resolve the smallest collapsed clumps 
of dark matter. This latter limitation can be severe when attempt- 
ing to predict the total annihilation radiation from a representative 
cosmological volume. For a standard supersymmetric neutralino, 
for example, the emission should be dominated by the smallest 
collapsed objects, w ith masses well below that of the Sun (e.g. 
iTavlor & S ilk 2003). Recent work has shown, however, that this 
problem is much less severe when predicting the observability of 
annihilation radiation from the Sola r System, which lie s just 8 kpc 
from the centre of the Milky Way (Springel et al. 200§j). Accord- 
ing to these authors, less than 3 percent of the dark matter within 
100 kpc of the Galactic Centre should be in small lumps; almost all 
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the rest should be in extended stream s of the kind discussed in this 
paper (see also lHelmi & Whitdl 19991) . They argue that the highest 
signal-to-noise for detecting the annihilation signal will be that of 
the smooth dark matter distribution in the inner few kiloparsecs of 
the Galaxy; small subhalos will be significantly more difficult to de- 
tect. Application of the techniques we have presented here should 
allow a rigorous evaluation of an important and previously unre- 
solved issue: whether the emission structure of this smooth compo- 
nent is significantly modified by caustic emission. In addition, they 
will allow an assessment of the expected morphology and observ- 
ability of outer caustics around external galaxies, most notably the 
Andromeda nebula. 

The authors thank Stephane Colombi, Craig Hogan, Roya Mo- 
hayaee and Volker Springel for helpful discussions. 
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